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Iterative algorithms with feedback are amongst the most powerful and versa- 
tile optimization methods for phase retrieval. Among these, the hybrid input- 
output algorithm has demonstrated practical solutions to giga-element non- 
linear phase retrieval problems, escaping local minima and producing images 
at resolutions beyond the capabilities of lens-based optical methods. Here the 
input-output iteration is improved by a lower dimensional subspace saddle- 
point optimization. 
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1. Introduction 

Phase retrieval is one of the toughest challenges in opti- 
mization, requiring the solution of large-scale, nonlinear, 
non-convex and non-smooth constrained problems. De- 
spite such challenge, efficient algorithms are being used in 
astronomical imaging, electron microscopy, lensless x-ray 
imaging (diffraction microscopy) and x-ray crystallogra- 
phy, substituting lenses and other optical elements in the 
image-forming process. 

In diffraction microscopy, photons scattered from an 
object (diffraction pattern) are recombined solving the gi- 
ant puzzle of placing millions of waves into a limited area. 
X-ray diffraction microscopy[l| has successfully been ap- 
plied to image objects as complex as biological cells 
quantum dots [H , nanocrystals [H and nanoscale aerogel 
structures Nanofabricated test objects were recon- 
structed computationally in 3D with several millions 10 
nm resolution elements [6] , other test patterns were cap- 
tured in the fastest flash image ever recorded at subop- 
tical resolution [7j. 

These experimental methods (see e.g. Q for a re- 
view) are being developed thanks to advances in op- 
timization techniques, primarily the introduction of a 
control-feedback method proposed by Fienup (Hybrid 
Input Output-HIO 0, Oil). The important theoreti- 
cal insight that these iterations may be viewed as pro- 
jections in Hilbert space [III E3] has allowed theoreti- 
cians to anal yze and improve on the basic HIO algo- 
rithm [H, H3, EE Ell • More recently Elser et al. 0|ll] 
connected the phase retrieval problem to other forms of 
"puzzles" , demonstrating performances of their Differ- 
ence Map algorithm 13 (a generalization of HIO) often 
superior to dedicated optimization pakages in problems 
as various as graph coloring, protein folding, sudoku and 



spin glass states. Rather than performing a local opti- 
mization of a constrained problem, the common theme of 
these algorithms is that they seek a solution to a different 
type of fixed point, the saddle-point of the difference be- 
tween antagonistic error metrics fl9| with respect to the 
feasible and unfeasible spaces defined by the constraints. 

Here each input-output iteration is improved by a 
lower dimensional subspace optimization of this saddle- 
point problem along the steepest descent-ascent direc- 
tions defined by the constraints. This lower dimensional 
optimization (performed here by Newton methods) is 
analogous to one dimensional line searches of gradient 
based methods, used to avoid overshooting and under- 
shooting in new search directions and providing faster 
and more reliable algorithms. 

The first sections introduce the phase retrieval problem 
and the saddle point optimization method, reformulating 
the HIO algorithm in terms of gradients and constraints 
(see [HI for further details). Sections [5] describe this 
lower dimensional optimization. Benchmarks performed 
on a simple simulated test pattern are described in the 
Sec. El 

2. Phase retrieval problem 

When we record the diffraction pattern intensity of light 
scattered by an object, the phase information is lost. 
Apart from normalization factors, an object with den- 
sity p(r), r being the coordinates in the object (or real) 
space, generates a diffraction pattern intensity equal to 
the modulus square of the Fourier Transform (FT) p(k): 



I{k) = \p(k)\ 2 



(1) 



'Current address: Lawrence Berkeley National Laboratory, 1 Cy- 
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where k represent the coordinate in the Fourier (or Re- 
ciprocal) space. In absence of constraints, any phase ip(k) 
can be applied to form our solution p = \f~Le %v . 

Phase retrieval consists in solving (Eq. |T|)) from the 
measured intensity values L(k) and some other prior 
knowledge (constraints). Diffraction microscopy solves 
the phase problem using the knowledge that the object 
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being imaged is isolated, it is assumed to be outside a 
region called support S: 

p(r) = 0, if r <£ S. (2) 

In practical experiments, where the object interacts by 
absorbing as well as refracting incident light, the prob- 
lem is generalized to the most difficult case of objects 
with complex "density" or index of refraction. In case of 
complex objects, this "support" region is the only con- 
straint and needs to be well defined, tightly wrapping the 
object. In many cases high contrast, sharp objects or il- 
luminating beam boundaries are sufficient to obtain such 
support of the object ab-initio @, H3| ■ 

A projection onto this set (P s ) involves setting to the 
components outside the support, while leaving the rest 
of the values unchanged: 

!p(r) if r G S, 
(3) 
otherwise. 

Its complementary projector can be expressed as P s = 

I Ps 

The projection to the nearest solution of (Eq. |T|)) in 
reciprocal space is obtained by setting the modulus to 
the measured one m(k) = y / I(k), and leaving the phase 
unchanged: 

P m p(k) = P m \p(k)\e lv{k) = Vl(k)e^ k \ (4) 

Singularities arise when p is close to 0, and a small change 
in its value will project on a distant point. Such pro- 
jector is a "diagonal" operator in Fourier space, acting 
element-by-element on each amplitude. When applied to 
real-space densities p(r), it becomes non-local, mixing 
every element with a forward T and inverse T~ x Fourier 
transform: 

P m = T- x P m T. (5) 
The Euclidean length of a vector p is defined as: 

\\p\\ 2 = ■ p = ^2\p(r)\ 2 = ^2\p(k)f- (6) 

r k 

If some noise a(k) is present, the sum should be weighted 
by w = -TT. The distance from the current point to the 
corresponding set \\Pp — p\\ is the basis for our error 
metrics: 

£ s O) = ||-P«P-P||, 

e m(p) = \\P m p-p\\, (7) 

or their normalized version e s ,m{p) = \\p m ^ P p\\ ■ 

The gradients of the squared error metrics can be ex- 
pressed in terms of projectors fiol |2~0] : 

Ve 2 m (p) = -2[P m -I]p (8) 
Ve 2 (p) = -2[P S -I]p, (9) 

Steps of — ^Ve 2 m bring the corresponding error metrics 
to 0. The solution, hopefully unique, is obtained when 
both error metrics are 0. 



3. Minimization in feasible space 

One popular algorithm [l(| H3| fits the data by minimiz- 
ing the error metric e m (p): 

mine 2 m {p), (10) 
p 

subject to Ps_p = 0, 

by enforcing the constraint and moving only in the fea- 
sible space p s — P s p. The problem is transformed into 
an unconstrained optimization with a reduced set of vari- 
ables p s : 

min£^(p s ). (11) 

Ps 

The steepest descent direction is projected in the feasible 
space: 

p (n+D _ p («) + Ap(«) 5 

ApW = -I Vs 4(>)), (12) 
= -P a [I-P m ]pW, 

where V s = P S V is the component of the gradient in the 
support. This algorithm is usually written as a projection 
algorithm: 

p (n+l) = PsPmp (n) . (13) 

by projecting back and forth between two sets, it con- 
verges to the local minimum. Such algorithm is com- 
monly referred to as Error Reduction (ER) in the phase 
retrieval community [10( . 

Notice that a step of — ■^7e 2 n {p) brings the error e 2 n {p) 
to 0. By projecting this step, setting to some of its 
components, we reduce the step length. Typically [l2l |23| 
the optimal step is longer than this step, we can move 
along this direction and minimize further the error met- 
ric. The simplest acceleration strategy, the steepest de- 
scent method, performs a line search of the local mini- 
mum in the steepest descent direction: 

min£^(p + (5Ap). (14) 



At a minimum any further movement in the direction of 
the current step increases the error metric; the gradient 
direction must be perpendicular to the current step. In 
the steepest descent case, where the step is proportional 
to the gradient, the current step and the next become 
orthogonal: 

^e 2 m ( P + SA Ps ) = (A Ps \2P s [I-P m }(p + SA Ps )) r , 
- (A Ps \[I-P m ](p + 6Ap s )) r , (15) 

where (x\y) r = (x' ■ y) and the projector P s is redun- 
dant and is removed. The line search algorithm can use 
e^j, and/or its derivative in (Eq. I|15p). This optimiza- 
tion should be performed in reciprocal space, where the 
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modulus projector is fast to compute (Eq. ([4}), while the 
support projection requires two Fourier transforms: 



P - TP S T- 



(16) 



but it needs to be computed only once to calculate Ap s . 

The steepest descent method is known to be ineffi- 
cient in the presence of long narrow valleys, where im- 
posing that successive steps be perpendicular causes the 
algorithm to zig-zag down the valley. This problem 
is solved by the non-linear conjugate gradient method 
0, HE HE US HI]- While error minimization methods 
converge to stable solutions quite rapidly however, the 
number of local minima is ofter too large for practical 
ab-initio phase retrieval applications and minimization 
methods are used only for final refinement. 

4. Saddle-point optimization 

The following algorithm is a reformulation of the HIO 
algorithm from a gradient/constraint perspective. We 
seek the saddlepoint of the error-metric difference C{p) = 



minmax£(p s + pg), 



(17) 



by moving in the steepest descent direction for p s 
(— P S V) and ascent direction (+jyv) for p L . For reasons 
discussed in appendix (|A"]) , we reduce the P s component 
by a relaxation parameter (3 S [0.5, 1]: 

V n > = {-P s +£Ps}|V£(>>). (18) 
The gradient of C (from Eqs. Jgjgj)): 

VC{ P ) = 2[P a -P m ]p, (19) 

is used in Eq. (JTHJ) to express the step and the new 
iteration point as: 



ApW = {P a [P m -T\-pP,P m }p 

[P s P m + Ps[I - PP m }}p 



(n) 
(«) 



(20) 



This iteration can be expressed in a more familiar form 
of the HIO algorithm 



expresi 

Urn 



(r) 



P m p {n) {r) 



if r £ S, 



(I - (3P m )p^ (r) otherwise. 



(21) 



Rather than setting to the object p(r) where it is known 
to be (r ^ S), this algorithm seeks a stable condition 
of a feedback system in which the nonlinear operator P m 
provides the feedback term P s P m p. From a fixed point 
whereby the feedback is but the constraint is violated 
{Ps_p ^ 0) it is often possible to obtain a solution by a 
simple projection P m p 13]. In fact P m p often satisfies 
the constraints better than the current iteration p, as the 
algorithm tries to escape a local minimum. 
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Fig. 1. Algorithms seek the intersection between two 
sets (data and constraints represented by two lines) us- 
ing different types of fixed points: ER projects back and 
forth between the two sets, moving toward the minimum 
of e m within the constraint (horizontal line). HIO seeks 
the saddle point (min s max,,) of (£), represented in back- 
ground in pseudo colormap. C is the difference of the 
square distances (e 2 ) between the current point and the 
two sets. The gradient (V£) is indicated by black ar- 
rows. In order to reach the saddle point, HIO spirals 
toward the solution by inverting the gradient component 
parallel to the constraint (horizontal direction) following 
a a descent-ascent direction ([— P s +Pg]V£ in light blue) 
toward the solution. Marker symbols are plotted every 
5 iterations. See also [T^ ] for a comparison with other 
algorithms. 



In the steepest descent method, optimization of the 
step length is obtained by increasing a multiplication fac- 
tor 6 until the current and next search directions become 
perpendicular to one another: 

(Ap\[-P s +pP s }VC(p + SAp)) r = 0. (22) 

A more robust strategy involves replacing the one di- 
mensional search with a two dimensional optimization of 
the saddle point: 



min max i/}(a, /3) , 

a f3 



ip{a, (3) = C(p + aA Ps + [3A P§ ) , 
A Ps = -|V s £(p);Ap, = |V,£(p); 



(23) 
(24) 



where both components (P s , Ps) of successive steps are 
perpendicular to one another: 



§| = (Ap a \V£(p + aA Pa +pAps)) r = 0, 
H = (ApJV£ (p + aA Pa + pA Pa )) r = 0. 



(25) 
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This two dimensional minmax problem needs to be fast 
to provide real acceleration, and will be discussed in the 
following section. 

5. Two dimensional subproblem 

The local saddle point (Eq. (f2"3")) ) requires two conditions 
to be met. The first order condition is that the solution 
is a stationary (or fixed) point, where the gradient of tp 
is (Eq. l[25])). We rewrite the condition in a compact 
form: 

\7 T ^(T) = (Ap\V p £(p + T T Ap) r = 0, (26) 
r = (2),Ap=(i£),V,= (£). (27) 

At the origin (r = 0), the gradient V T ip is negative in the 
P s subspace and positive in the subspace, decreasing 
(increasing) C in the steepest descent (ascent) directions 
in the two orthogonal spaces: 

v^(o) = (28) 

= „ /H|Ap 3 || 2 \ _ i /-||V s £(p)|| 2 \ 
A \+\\&P±\\ 2 ) 2 \+\\V ± C(p)\\ 2 ) ■ 

The minimal residual method finds a stationary point by 
minimizing the norm of the gradient: 

inin$(r), $= \\\V^j(t)\\ 2 , (29) 

transforming the saddle point problem in a minimiza- 
tion problem, and providing the metric $ to monitor 
progress. However this method can move to other sta- 
tionary points. 

Second order conditions (min-max) require the Hessian 
H of -0 (the Jacobian of (Eq. (f2"6"|) ) to be symmetric and 
indefinite (neither positive nor negative definite): 

^(sa)*)^: (30) 

This Hessian is computed analytically (see appendix ([Bj , 
it is small (2x2), and can be used to compute the Newton 
step: 

At = -W _1 V t V>. (31) 

However, the Hessian precise value is not necessary and 
requires an effort that could be avoided by other meth- 
ods. 

The normalized steepest descent-ascent (SDA) direc- 
tion can be expressed in terms of a Newton step using an 
approximate diagonal Hessian whose elements are equal 
to -V</>(0): 

-O.SDA _ 9 / l|Ap s || 2 \ ,„ 9 \ 

The Hessians Ti satisfies condition (Eq. (|30p). ensuring 
that At is less then 90° from the direction of the saddle. 
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Fig. 2. Pseudocolor and contour maps of the lower 
dimensional function ip(a, 0) are depicted in the back- 
ground. This ip was computed from one iteration of the 
test described in Sec. (JSJ- The saddle, typical of this 
test, is wide in the horizontal direction and narrow in 
the vertical, showing the importance of the relaxation 
parameter (3 < 1 to reduce (precondition) the vertical 
move. However, the condition number increases near the 
saddle-point as it becomes even narrower in the verti- 
cal direction. Often the dependence of ip with respect 
to the fitting parameter (3 (vertical axis) resembled a v 
centered near the solution, and successive iterations us- 
ing just a preconditioner jumped up and down near the 
saddle-point. Iterations using the SRI quasi-Newton up- 
date (starting from — 0, TL^ — 7i HI °) shown here 
converge rapidly to the solution. 

Starting from = 0, the first iteration gives a unit 
step = — H~ 1 Vip(0) = (I). A preconditioner can 
be used to reduce the feedback: 

n mo = (Sx/°,-)w SDA , (33) 

providing the HIO step at the first iteration, = 
At' ' = (jj)- This approximate Hessian can be used as 
a starting guess, which is often all it is needed to achieve 
fast convergence. 

We can perform a line search using the preconditioner 

^At|W _1 V x ^(t + SAt)^ =0. (34) 

However the Hessian of HT 1 ^ is antisymmetric, the al- 
gorithm is unstable and could spiral away from the solu- 
tion. The bi-conjugate gradient method applies to sym- 
metric indefinite Hessians and monitors progress of the 
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algorithm. Conjugate directions At replace the steepest 
descent direction in the line search (with y(°' = 0): 



Ar («+i) 

7 («) 



Ar (»+1) + 7 (n) Ar (n) ^ (35) 

(AT< n + 1 '|H(- 1 )(V^(r<"+ 1 ))-VV(r<'«)))) 
(At-(")|'H- 1 AVV(i-<"))) ' 



A better option is to use a quasi-Newton update of the 
Hessian or its inverse (a secant method in higher dimen- 
sions) based on the new gradient values. The Symmetric 
Rank 1 (SRI) method can be applied to indefinite prob- 
lems 1291: 




V 

ah- 1 
n- 1 



Fig. 3. Test figure used for benchmarking (total size: 
256 2 , object size: 128 2 pixels. The support was slightly 
larger than the object: 129 2 pixels. 



(At 

n- 1 - 



H-iy) T -y 
AH' 1 . 



(36) 



t + At and go back to 4. 
p+T T Ap. If e m is small exit, otherwise 



Second order conditions (Eq. ([20])) can be imposed to 
the Hessian, by flipping the sign or setting to the values 
that violate them. $ can be used to monitor progress, as 
long as we are in the neighborhood of the solution and 
the Hessian satisfies second order conditions (Eq. (I30D V 
It was found that the Hessian and step size parameters 
where fairly constant for each 2D optimization, therefore 
the first guess for r and TL was obtained from the average 
of the previous 5 of such optimizations. With such initial 
guess, 3 SRI iterations of the lower dimensional problem 
where often sufficient to reduce $ below a threshold of 
0.01 $(0) in the tests described below. 

In summary, an efficient algorithm is obtained by 
a combination of HIO/quasi-Newton/Newton methods 
with a trust region |At| < r: 



9. update r 

10. update p - 
go back to 1. 

The trust region is used to obtain a more robust algo- 
rithm, is reset at each inner loop, it increases if things 
are going well, decreases if the iteration is having trou- 
ble, but it is kept between (r m ; n , r max ), typically (0.5, 3). 
We can keep track of t, Vt computed, and restart the 
algorithm once in a while from the root of the 2D linear 
fit of V>(t). 

We can easily extend this algorithm to two succes- 
sive steepest descent-ascent directions, by performing a 
higher dimensional (4D) saddle-point optimization: 



min max C.L 



Jn-i 



l))T A/9 (n+l) +(r (n) ) T Ap 



(n) 



1. calculate step Ap = — | ( 
region radius r = r max . 



-V S C h 



and set trust 



This 4D optimization is performed using the same 
Newton/quasi-Newton trust-region optimization as in 
the 2D case. The first step (r°) T Ap° is obtained solving 
the 2D minmax problem, and the following ones will be 
perpendicular to the last 2 iterations. 



2. if the iteration number is < 5, use HIO as first 6. Performance tests 



guess: TL = H 



HIO 



-(i) 



(I,/?)- 



3. otherwise average 5 previous optimized step sizes 
t, and Hessians H, and use the average as initial 
guess. 

4. calculate gradient Vt/>(t). If small, exit loop (go 

tonoj). 

5. compute Newton step using approximate Hessian: 
At = — 7i~ 1 V , 0, enforce trust region |At| < r. 

6. update inverse Hessian with SRI method (Eq. [55]) . 

7. if the Hessian error 1 1 At — 7Y _1 i/| | 2 is too large, 
calculate the true Hessian, perform a line search, 
decrease trust region radius r. 

8. force Hessian to satisfy second order conditions 
(Eq. [317]) . by changing the sign of the values that 
violate conditions. 



The tests presented here were combined with other al- 
gorithms in [l9l |. In such comparison between algo- 
rithms, the original HIO algorithm performed better than 
other algorithms described in the literature. Fig. [3] was 
used to simulate a diffraction pattern, and several phase 
retrieval tests were performed using different random 
starts. When applying a nonnegativity constraint HIO 
always converged within a few hundred iterations. The 
algorithms were therefore tested using a more difficult 
problem, nonnegativity and reality constraints were re- 
moved, allowing the reconstructed image to be complex, 
adding many degrees of freedom within the constraint. 

When the error metric e m fell below a threshold it 
was considered a successful reconstruction. The thresh- 
old 10 -4 chosen was enough to obtain visually good re- 
constructions. 

Fig. 0] shows the relative performance of the various 
algorithms. By adding 2D or 4D optimization, the al- 
gorithm converged more reliably and in less iterations 
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100 



— ^J=o -a 




tA-HIO 

-D-S02D 

-O-S04D 





5000 10000 

# iterations 



Fig. 4. Succesfull phase retrievals, starting from random 
phases, as a function of iterations. When the normalized 
error metric e m falls below a threshold (10 -4 ) it is con- 
sidered a succesfull reconstruction. The plot represents 
the cumulative sum of the number of succesfull recon- 
structions vs the number of iterations required. Values 
for 50% (dotted line) and 100% success are listed in Tab. 
[TJ Comparison with other algorithms is shown in (l9j . 

Table 1. Benchmark of various algorithms (250 trials) 
summcrizing the results in Fig. |4] 



Algorithm 


No. of iterations for 


success rate 




50% success 


100% success 


(total) 


HIO 


2790 


> 10000 


82% 


S02D 


656 


5259 


100% 


S04D 


605 


2999 


100% 



to a solution (Table [TJ). The 2D (and 4D) optimization, 
written in an upper-level language (matlab) increased the 
computational time of each iteration by a factor of 3 (and 
4) and required the storage 2 (and 4) additional matrices 
compared to HIO. HIO itself requires 2 matrices in addi- 
tion to the data and constraints. A c-code version of the 
algorithm developed by F. Maia from U. Uppsala per- 
formed slightly faster, and could be further optimized by 
reducing some of the redundant computations involved. 

This lower dimensional optimization employs matrix 
cross products, computing a number or floating point 
operations proportional to the number of elements of the 
images (and can be implemented on parallel systems). 
The Fourier transforms employed to calculate the steps 
in the higher dimensional problem will dominate the com- 
putational burden for larger matrices. 

7. Conclusions 

The hybrid input-output method, usually described as 
a projection algorithm, is a remarkable method to solve 
nonconvex phase problems. When written as a saddle 
point problem, HIO speed and reliability can be improved 



by applying Newton methods to explore lower dimen- 
sional search directions. This approach differs from other 
nonlinear optimization algorithms that try to satisfy con- 
straints using various forms of barriers or trust regions, 
often requiring stochastic methods to climb out of local 
minima present in nonconvex problems. 

The saddle-point algorithm, while following a path in- 
dicated by a gradient, does not seek a minimum and does 
not stop at local minima. Although stagnation occurs, 
it appears that the area of convergence to the global so- 
lution is much larger compared to simple minimization 
methods. Problems arise at the locations of intensity 
values where phase singularities occur, causing optimiza- 
tion problems to be nonsmooth. Various forms of pre- 
conditioning low intensity values could be applied to this 
problem [30L l3l| . here the nonsmooth behavior of the 
problem is addressed using a lower dimensional optimiza- 
tion to stabilize the iterations, providing a more reliable 
algorithm. 

In this paper, the saddle-point optimization was per- 
formed for a simple equality constraint P L p = 0. Such 
linear constraint allows rapid calculation of the gradients 
involved in the lower dimensional optimization. Linear 
approximation or more advanced methods could be ap- 
plied to other nonlinear constraints such as thresholds, 
histograms, atomicity and object connectivity, extending 
this approach to a larger class of problems. This saddle- 
point optimization formalism can easily be generalized to 
other problems of conflicting requirements where gradi- 
ents or projections can be computed. 

Appendix A: Relaxation parameter and phase 
singularities 

The large dimensional minmax problem (Eq. I17|) can be 
expressed in a system of two parts: 

! min p s £ li(P) = mhl Ps HI 7 - P m]p\\ 2 
(Al) 
min p, £ 2 s (p) - £ 2 m (p) = min Pi 2(P ro p|p) r + c 

The upper optimization is similar to the problem treated 
in Section^ converging to a local minimum with a simple 
projected gradient method. The lower function, however, 
can be discontinuous in the presence of zeros (p s = 0) in 
Fourier space: 

(P m p\p) = ]T Vl\p s + Ps} (A2) 

which is a non-smooth v-shaped function of p s for p s — 
0, \fl > 0, and simple gradient methods oscillate around 
minima. The projected gradient step can be overesti- 
mated and requires the relaxation parameter (3. Zeros in 
Fourier space are candidates (necessary but not sufficient 
condition) for the location of phase vortices, phase dis- 
continuities, which are known to cause stagnation [33|. 
Analytical [32j , statistical [1, HI, [34[ , and deterministic 
[3fl l3l| methods have been proposed to overcome such 
singularities. 
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Appendix B: Two dimensional gradient and Hes- 
sian 

The function C in reciprocal space can be expressed as: 



C(p a +ps) = - P m ](ps + Ps) -\\p.\\ ( B1 ) 

= \p> I 2 - 2v7|p s + psl + V7. 

and the two components of the gradients: 

vr _ ( p s vc \ _ „ ( P s [I-P m ](Ps+Ps) \ (T >o\ 

VL ~ [Psvcj-ty -p,p m (p e+Pj! ) ) \ al > 

The corresponding steps Ap s = — |V S £, Ap ± = +^\7 ± £. 

The function ip(a,/3) can be calculated in recipro- 
cal space, provided that the components p Sj£ , Ap s _s are 
known: 

i>(a,(3) = \\[I-P m ](p + aAp s +f3Ap s ) 

- \\Ps +0Aps\\ 2 
k 

l 2 

- | p £ + [3Aps] 

= ^2 1 + \p s + aAp s | 2 
fe 

-2\/7|p + aAp s + /3Ap s | 



(B3) 



The gradient components (writing p T = p 
p + aAp s + fiAps) are: 

VtV = (Ap|V£(p T )) r , 

= 2(Ap|[P s -P m ]p T ) r , 

^ = 2(Ap 3 \[I-P m ]p T ) r , 

% = 2(Ap s _\[-P m ]p T ) r , 

Using common derivative rules: 

d 



r T Ap 



(B4) 



dx 



da I 



a Ax | 

9 xjaAg 
da \x+aAx\ 



R^Ax^ (x+aAx)) 
\x-\-aAx\ 
Ax x-\-aAx 



(B5) 

(B6) 
(B7) 



£\x + aAx\ 2 
x + aAx\ 



da 

d 2 
da' 2 



|x+aAa;| |:c-t-aAx| 

i R( Ax' ! (x+aAx)) 
\x+aAx\ 2 ' 

23? (AxHx + aAx)) , 

\Ax\ 2 _ U(AxHx+aAx)) 2 
\x+aAx\ \x+aAx\ 3 ' V 13y / 



^- l x + aAx\ 2 



\Ax\ 2 , 



and 



d 2 \x+aAx+pAy\ }R(Ax' 1 Ay) 

dad/3 ~ \x+aAx+f3Ay\ 

SR(Axt (x+aAx+PAy)) 
~ \x+aAx+[)Ay\ 

SR(Aj/t (x+aAx+pAy)) 
\x+aAx+f3Ay\ 2 ' 



(B8) 
, (B9) 
(BIO) 

(Bll) 



we can calculate the analytic expression for the Hessian 
using P m p = ||| \fl. Starting from the simplest compo- 
nent: 



a 2 V 



-2<ApV^^ 



9/3 



(B12) 



, |Apj 2 V7 , SR(Apj,pv) 2 y? 

_ V- |ApJ 2 y7 v7 ^(Apj,pV)(Apj,pV+Ap,pt) 

->(*y-*(l-4w)IAA} r 



9jj/> 
da 2 



2(Ap. 



- I 0[J-P m 1p T 



9a 



(B13) 



= 2^|Ap s | 2 -^ + ^S^ 



= 2(Ap. 
The cross terms: 



I 1 ~ W7\ ( 



Apf ' 



Ap, 



dpda 



2(Ap s 



= 2^Ap s 
= 2/Ap, 



9[/-P m ]p T 

9/3 

2|p T | 

2WI ^ 



(B14) 



Apf |pVI 2 
# |Ap,| 2 



Ap t 
A~ Pt 
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